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An explicit, practical procedure is suggested for transforming from the laboratory variables density (p) and 
temperature (T) into the parametric variables r and 0, which occur in various scaled representations of equations 
of state and of transport properties of fluids near critical points. A reasonably efficient and versatile computer 
program illustrating this procedure is provided. With this program, the parametric equations of state which occur 
in several formulations of simple, extended, and/or revised scaling are as easy to use as any other equation of state 
for which T and p are the independent variables. 
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1 . Introduction 

In recent years various scaling equations of state have 
been developed to provide theoretically based, highly accu- 
rate representations of thermodynamic data for pure fluids 
and fluid mixtures near their critical points. These equations 
of state are all consistent with the hypothesis that very near 
the critical point the thermodynamic potential is a homoge- 
neous function of its natural variables. These equations of 
state have relatively few parameters which must be adjusted 
to fit data for any specific fluid. This is a practical advantage 
for those concerned with taking or correlating data. On the 
other hand, many scaling equations of state are written in 
terms of parametric variables which appear in nonlinear 
equations linking together such physically relevant variables 
as the pressure (P), the density (p), and the temperature (7"). 
The occurrence of the parametric variables is a major 
practical disadvantage to any nonexpert who might otherwise 
be interested in using a scaling equation of state to deal with 
a limited problem. It is, of course, true that these parametric 
variables can be eliminated algebraically from the equation 
of state along special paths in thermodynamic space such as 
the vapor pressure curve. Then, along the special paths, the 
scaling equations reduce to simple power laws with noninte- 
gral exponents relating directly measured quantities to each 
other. For states off these special paths, a user of a scaling 
equation of state must eliminate the parametric variables 
numerically. The purpose of the present paper is to provide 
an explicit example of one way in which this can be done. If 
our method is followed, a user need not repeat the program- 
ming effort required to eliminate the parametric variables in 
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those situations in which p and T are independent variables. 
In situtations in which it is convenient to designate other 
physical variables as independent, the general method out- 
lined here may still be followed. 

In this work we discuss only those scaling equations of 
state which use parametric variables. There are several 
published examples of scaling equations of state which do 
not use parametric variables [1,2]. 3 Unfortunately, these 
non-parametric equations of state cannot be integrated in 
closed form to obtain an explicit expression for the thermo- 
dynamic potential. Thus, frequently used expressions for the 
entropy and for the specific heats consistent with non- 
parametric equations of state contain integrals which must be 
evaluated numerically. In our experience, the lack of an 
explicit thermodynamic potential has been a serious handi- 
cap in attempting to use these non-parametric equations for 
constructing thermodynamic models for real fluids. Accord- 
ingly, we now choose to deal with the numerical problems 
arising from the use of parametric variables rather than the 
numerical problems associated with integral representations 
for the potential. 

2. Simple Scaling with the Linear Model 

We will consider first the "linear model" parametric 
equation of state proposed by Schofield [3]. This equation 
was first proposed on phenomenological grounds and was 
applied with success to the insulating ferromagnet CrBr 3 [4]. 
It has since received some theoretical justification [5] and it 
has been widely used to correlate equation of state and 
transport property data in various fluid systems [6-9]. Most 
important from our present point of view, the equations to be 
solved in using the linear model parametric equation of state 
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may be easily generalized to handle the more complicated 
cases of extended and/or revised scaling. 

To use the linear model parametric equation of state in 
situations where the temperature (T) and the density (p) are 
the independent variables, the parametric variables r and 8 
must be obtained by solving the following two equations: 



T = T c + 7X1 - b 2 2 ) 

P — Pc + Pc mrfi Q 



(1) 

(2) 



Here T c and p c are critical constants. For the simple fluids 
studied by Levelt Sengers et al., [7] misa positive constant 
in the range 1.5 < m < 2.1; b 2 is a constant falling in the 
range 1.1 < b 2 < 1.6, and the exponent /3 usually falls in 
the range 0.31 < (3 < 0.38. The transformation represented 
by eqs (1) and (2) for the fluid xenon is sketched in figure 1. 
To make this sketch, values of T c , p c , m, and b 2 were taken 
from reference [7]. Once values of r and are obtained by 
solving (1) and (2), other quantities such as the pressure 
(P) and the constant volume specific heat (C v ) may be 
calculated as functions of p and T by simply substituting r 
and 6 into algebraic equations which can be solved explicitly 
for the desired quantities [7,10]. There are many ways to 
solve these equations numerically. The method described 
here is reasonably efficient and reliable. 

To solve (1) and (2) we first introduce the dimensionless 
variables A7 7 * = (T - T C )/T C9 and Ap* = (p - p c )/p c to 
obtain 



o = Ar* - r(i - b 2 e 2 ) 

= Ap* - rnr^O 
We then eliminate r from (1) and (2) to obtain 



= - Ap* - 1 Ar* j 

m 



-e + be 1 1 - b 2 e 2 1 



(3) 
(4) 

(5) 



Note that absolute value signs had to be introduced because 
we have chosen to raise to powers the separate quantities 
AT* and (1 — b 2 6^) rather than their quotient which always 
remains positive. By introducing two new symbols, Z and C 



z = be 



C = -Ap* - I A7 7 * I 
m 



(6) 



-a 



we obtain a compactly written transcendental equation in one 
variable 



f(z) = = c + z\i - z 2 \-e 



(7) 



Numerical approximations to its roots may be found effi- 
ciently by Newton's method. This is particularly easy to 



implement because the derivative, /'(Z), may be calculated 
analytically. In Newton's method, the (n + l)'th estimate for 
Z is found from the n'th estimate by the rule 



Zn+i ~ Z n 



f(Z) 
f(Z) 



In the present case this rule takes the form 

(l-Zl)(Z n + C n \l-Zl\ B ) 



^rtA-} ^n 



(1 - Zl) + 2(3Zl 



(8) 



(9) 



Equation (9) is quite satisfactory for computation if suitable 
precautions are taken to handle the singularity at Z = 1. 
Each iteration requires only a single exponentiation. A 
reasonable initial estimate for Z is obtained by looking at the 
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FIGURE 1 . Sketch of the transformation from p and T to r and 6 in the 
linear model (eqs (1) and (2)) with parameters appropriate for xenon from 
reference [7]. 



physically important special curves, namely: the critical 
isochore above T c for which Z = 0, the critical isotherm for 
which Z = ± 1 and finally or coexistence curve on which Z 
= + b in the liquid phase and Z = — b in the vapor phase. 
An estimate which has the appropriate limits near the special 
curves is: 
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Z = 1 



(1 - b)AT* 



1 - b 2 
Z = (1 + Ar* 



Ap" 



1/0 



6Ap* 



for AT* <0 



-P for AT* > 



(10) 



We have applied this approach to the solution of eqs (3) and 
(4) with the double precision Fortran program listed in 
appendix I. For all of the physically reasonable values of 
input parameters we have tried, at most five iterations 
were needed. The resulting values of r and 8 satisfy eqs (3) 
and (4) to 12 or more significant figures. (When |Ap*| and/or 
\AT*\ are smaller than 10~ 12 , an "unphysical" situation, the 
sample program yields fewer significant figures). If a sub- 
stantial amount of computation were to be done with a single 
value of j8, a program could be written in which the time- 
consuming step of raising a floating point number to a 
floating point power occurs less frequently than the 2-7 
times it occurs with the present program. 

In concluding this section we note that solutions to (3) and 
(4) in the range 1 < 8 2 < l/[/> 2 (l - 2/3) ] may be interpreted 
as either two phase states (for which 6 should be set to ± 1) 
or metastable extensions of the single phase isotherms [11]. 
Very little is known about the accuracy of the linear model 
when it is used in the metastable region. 

3. Additional Examples Using Simple Scaling 

Other familiar parametric versions of simple scaling can 
be handled in the same way as the "linear model". For 
example, the "cubic model" introduced by Ho and Litster 
[12] uses two parametric variables which may be considered 
to be defined by the equations: 



= AT* - r(l - b 2 2 ) 
= Ap* - mr p 0(l + c 2 2 



(ID 
(12) 



Again we have used r and 8 as the parametric variables while 
m, b, (S, and c 2 are constants. If the same substitutions are 
made in (11) and (12) as those used above in connection with 
the linear model, one may obtain: 



f(Z) = = C + Z(l +gz 2 )| 



I 



(13) 



In practice, c 2 Z/b 2 is much less than unity so that (13) is 
very similar to (7). Initial estimates for the solution to (13) 
may be obtained from eq (10), just as for the "linear model"; 
then the solution by Newton's method proceeds exactly as 
outlined above except that the "cubic model" function for 
f(Z) is used instead of the linear model function. 



A parametric scaling equation of state was introduced by 
Wilcox and Estler [13] to deal with data from diffraction 
experiments in fluids near critical points. Their scheme has 
been used by Estler et al. [14] and by Hocken and Moldover 
[15]. For a given fluid state, the parametric variables in the 
Wilcox-Estler scheme have numerical values which are quite 
different from those which appear in the "linear model" or 
the "cubic model". One may consider the variables to be 
defined by the equations: 



= \Ap*\-[Y R(l-d/8 o ry 



= A7'* 



I + 



R0 



j8 i - e/e, 



y i - e/e e _ 



(14) 



(15) 



Here, following Wilcox and Estler, the parametric variables 
are "/?" and "0" while F , O , 8 X , /3 and y are all constants, 
and A is a combination of constants: A = 1 — 8J8 X . (Note: 
Because the absolute value of Ap* appears in (14), the 
inverse transformation from R and 8 to Ap* and AT* is not 
unique). Once again one of the parametric variables, say /?, 
may be eliminated from (14) and (15) to yield an equation in 
one variable which may be solved by Newton's method. 

On occasion, T and p are not the most convenient 
independent variables. For example, the analysis of flow 
calorimetry data [16] requires the computation of the en- 
thalpy as a function of T and P. Each of the simple scaling 
equations of state mentioned above will lead to two equations 
relating P and T to parametric variables. The numerical 
problem of solving these equations for r and 8 is quite similar 
to the problem we have discussed above, namely of solving 
the equations in p and T for r and 8. In particular, when P 
and T are known it is still possible to eliminate one of the 
parametric variables from the two equations algebraically. 
Thus a single equation in one unknown remains which may 
be efficiently solved with Newton's method. 

4. Extended and/or Revised Scaling 

In more complicated versions of scaling, the conversion 
from the laboratory variables p and T to the parametric 
variables r and 8 requires the numerical solution of two 
simultaneous equations. One parameter can no longer be 
eliminated algebraically. An efficient and versatile approach 
is to use simple scaling to obtain a first approximation to the 
values of the parametric variables. The approximation may 
then be improved by using Newton's method for the solution 
of simultaneous equations. 

The same first approximation (simple scaling) may be used 
with a variety of more complicated models. Where the 
derivatives required for the use of Newton's method are 
computed numerically, it becomes quite easy to change from 
one model to another. 
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We illustrate this approach by considering the form of 
extended and revised scaling used by Balfour et al. [10] to 
study the equation of state of steam near its critical point. In 
this scheme the defining equations may be written in a form 
analogous to eqs (3) and (4) above: 



C/(r,0) - = A7 1 * - r(l - b 2 2 ) 



(16) 



1 — qar 
V(r,6) = = Ap$ - rnr^S 



1 + — + qar l ~ 
m 






' S (6)/6 



(17) 



+ qak] 



r l-a-/3+A 



Sm/imB)} 



Here, the variable Ap$ is the scaled density deviation from 
the "rectilinear diameter" (In the notation of Balfour et al. 
Ap$ = 4?i A7 1 * + Ap*). This contrasts with Ap* used 
above which is the scaled density deviation from the critical 
density. In (16) and (17) the exponents a and A are about 
0.1 and 0.5, respectively. The constants q and kjm are 
small compared with unity, and S (6) and Si(0) are rational 
polynomial functions of 6. 

To consider the transformations (16) and (17), we recall 
that both revised and extended scaling may be cast in the 
form of an expansion about the critical point. Thus the small 
"revision" coefficient q is a measure of the lowest order 
departure of the isotherms from antisymmetry and the small 
"extension" coefficient k 1 /m is a measure of the strength of 
the second most singular term in the free energy expansion 
about the critical point. We expect that if the expansion is to 
make sense, the terms with these coefficients should be 
small compared with the leading terms. Thus, for physically 
meaningful values of Ap$ and A7 7 *, good initial estimates 
for the parametric variables r and 6 in revised and/or 
extended scaling can be found by setting k l /m and q equal to 
zero. This, of course, reduces eqs (10) and (11) to eqs (3) 
and (4) which we have already solved in dealing with simple 
scaling. 

Then, the initial estimates for r and 6 may be improved 
efficiently by using Newton's method for the solution of 
equations in two variables. The (n + l)'th estimates for r 
and 6 are obtained from the ra'th estimate by the rules [17] 



r n+i ~ r n 



Qn+l - On 



,W , , dU 
U(r n ,6 n ) — - V(r n ,O n ) — ID (18) 



V(r n M y- ~ U(r n ,6 n ) j- /D 
or dr 



BUdV dUdV 

D = - — — (20) 

dr dS dO dr 



The partial derivatives in eqs (18), (19), and (20) are 
understood to be evaluated near (r n , 6 n ). In practice we have 
computed approximations for them from divided differences. 
Thus, for each iteration, one must evaluate U and V at three 
sets of coordinates: (r,0), (r 4- 8,9), and, (r,0 + 8). The 
small difference, 6, used for estimating the derivatives is 
chosen to be 10 -8 in the sample program in the appendix. 
This value is satisfactory for all the physically reasonable 
values of AT* and ApJ we have examined. The singularity 
at A7 1 * = ApJ = does not cause numerical problems 
because, near this point, simple scaling gives quite accurate 
values of r and 0, and the iteration of eqs (18-20) is not 
used. 

If one wished to alter the scheme of Balfour et al., say by 
adding an additional extension to scaling term with a new 
coefficient, /c 2 , the only change needed in this computational 
scheme woud be the addition of that term to the functions 
U(r,6) andV(r,d). 

It is economical to evaluate U and V in the same 
subprogram to minimize the number of exponentiations 
required in each iteration. 
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6. Appendix I 

The subroutine RTHETA solves eqs (3) and (4) of the text 
for the parametric variables r and 6. The argument list of the 
subroutine eorresponds to the symbols in those equations in 
the order: r, 0, Ap*, A7 1 *, /3, m, b 2 . This version treats 
thermodynamie states for whieh 1 < | 6 | < 1.00234 as 
though such states were metastable. (The number 1.00234 is 
arbitrarily chosen. It must fall within the range indicated in 
the text). States for which 1.00234 < \ 0\ are treated as two 
phase states. In other words, 6 is assigned the value + 1 or 
— 1 which is appropriate for an equilibrium Liquid or vapor 
state at coexistence conditions. 



SUBROUTINE RTHETA (R . THETA . RHO . TEE • BETA, EM. BESQ ) 
IMPLICIT REAL*8 (A-HtO-Z) 
450 FORMAT(2X^ f CATA ERROR IN RTHETA* 9 7(1X. E12.6) ) 
IPC EM.LE. 0.D0 .OR. BESQ • LE • 1*00 ) GO tq 600 
ABSRHC = DABS(RHO) 

IF C ABSRHO .LT. l.D-12 ) GO TO 600 
BEE a DSQRTC BESQ ) 

IF CDABS(TSE).LT«1«0D-12 I GO TO 495 
IF ( TEE .LT. O.DO ) Z = 1.00 - ( 1. DO-BEE ) 
1 * TFE / (1.D0-BESQ) * ( EM/ABSRHO ) **< 1 .DO /BETA ) 
IFCTFE.GT. O.DO) Z -( 1 •D0^TEE*< FM/BEE/A8SRHQ )* * ( 1 .DO/ BETA) )**-BETA 

IF ( Z .GT. 1.0J234D0*BEE ) GO TO 496 
C a -RHO*Bf£/EM/CABS<T£E)**BFTA 
Z a DSIGNCZ.RHO) 
100 DO 500 N = 1» 10 
12 a Z*Z 
Z3 = 1.00 - Z2 

DZ a Z3*(Z ♦ C* DABSC Z3)**BETA)/< Z3 + 2 • D0*8ETA* Z2 > 
Z a Z - DZ 

IF ( DABSC DZ/Z I .LT. l.D-12 ) GO TO 498 
500 CONTINUE 
601 WRITEC6.450) R • THETA .RHO • TEE. BETA , EM , BESQ 

RETURN 
498 THETA a z / BEE 

R = TEE / ( 1.D0 - Z*Z ) 
RETURN 

495 THETA = DSIGNC 1. DO. RHO ) / BEE 

p s ( RHO / ( EM * THETA ) ) ** (1. DO/BETA) 

RETURN 

496 THETA = OS I GN ( 1 .DO. RHO ) 

R = TEE / ( l.ODO - BESQ ) 

RETURN 
600 IF ( DABS(TEE) .LT. l.D-12 ) GO TC 601 
IF( TEE .LT. O.DO ) GO TO 496 
THETA a O.DO 
R a TEE 

RETURN 
END 
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7. Appendix II 

This subroutine, SOLVE, solves equations similar to (16) 
and (17) of the text for the parametric variables r and 0, The 
argument list is identical with that in appendix I. This 
subroutine assumes that the subroutine RTHETA of appendix 
I is available and that another subprogram, UV, is available 
to evaluate the functions U(r,6) and V(r,6) defined by 
equations such as (16) and (17). The UV subprogram's 
argument list corresponds to the symbols r, 0, U(r,0), V(r,6), 
Ap*, A7 7 *, as they are used in eqs (16) and (17). Statement 
10 of SOLVE insures that 6 will be assigned a value of + 1 or 
— 1 if the values of Ap* and A7 1 * provided are such as to lie 
within the coexistence curve. Thus this version of SOLVE 
does not assume the existence of metastable states. 



SUBROUTINE SOLVEC R .THETA • RHO, TEE. BETA* EM ♦ BESQ) 
IMPLICIT REAL*8 (A-H.O-Z) 

IFCTEE.GE. Q.DO) GO TO 55 

THETA = DSIGNC 1 .D0.RHO) 

R = TEE / C 1 .CO-BESQ) 

CALL UV{*,THETA.U.V.RHO.TEE ) 
10 IF( (V*RHO) ,LE. 0-00 ) GO TO 999 
55 CALL RTHE7A<R. THETA. RHO, TEE. BETA. EM. BESQ) 

DO = l.D-8 
65 DO 75 I = lt6 

CALL UV(R. THETA. U. V.RHO»TEE> 

IF(0ABS(V) .LT. l.D-5 .AND. DABS(U) .LT. l.D-5) GO TO 999 

CALL UV C R+CC .THETA 9 UR t VR .RHO. TEF > 

CALL UV( R, THETA +DD.UT.VT. RHO. TEE) 

DUDR a UR - U 

DUOT sUT- U 

DVDR = VR - V 

DVDT a VT - V 

D as OUCR*DVDT - DUCT*DVDR 

R = 9 - DD*(U*DVDT - V*DUDT)/D 

THETA = THETA - CD*|DUDR*V - DVDR *U ) /D 
75 CONTINUE 
998 WRITEC6.671 RHO. TEE. R. THETA. BETA. EM. BESQ 
67 F0RMATC2X. • SOLVE .DOES NOT CCNVERGE f . 7F15.8) 
999 CONTINUE 

RETURN 

END 
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